%write annual harmonics of ocean velocity
%as input to alexei's program
%3dec07 as pointed out by stefan, a testrun is designed
% first keep all imag as zeros and rund x3dg with > 4000days
% second keep all reals zero and run x3dg with > 4000 days
load c:\manoj\ocean\annual_ecco_0006;
load c:\manoj\ocean\ECCO_Cord;
[X1,Y1] = meshgrid(V_long,U_lat);

IUr = griddata(U_long,U_lat,real(data_iu),X1,Y1,'linear');
IUi = griddata(U_long,U_lat,imag(data_iu),X1,Y1,'linear');
IVr = griddata(V_long,V_lat,real(data_iv),X1,Y1,'linear');
IVi = griddata(V_long,V_lat,imag(data_iv),X1,Y1,'linear');
IUr(:,360) = IUr(:,1);
IUi(:,360) = IUi(:,1);
IVi(160,:) = 0;
IVr(160,:) = 0;

IUr_n = [zeros([10,360]);IUr;zeros([10,360])];
IUi_n = [zeros([10,360]);IUi;zeros([10,360])];
IVr_n = [zeros([10,360]);IVr;zeros([10,360])];
IVi_n = [zeros([10,360]);IVi;zeros([10,360])];

data_iu = complex(IUr_n,IUi_n);
data_iv = complex(IVr_n,IVi_n);

load C:\Manoj\projects\tsunami\Bz;
Bz = rot90(Bz);
data_ur = flipud(real(data_iu))*0;
data_ui = flipud(imag(data_iu));
data_vr = flipud(real(data_iv))*0;
data_vi = flipud(imag(data_iv));


prd = [365]; % The mean flow has very high period -> 4000 days !

% fp = fopen('d:\manoj\temp\a.project','wt');

  
CUr = -data_vr.*Bz.*1e-9.*3.2;
CUi = -data_vi.*Bz.*1e-9.*3.2;
CVr = -data_ur.*Bz.*1e-9.*3.2;
CVi = -data_ui.*Bz.*1e-9.*3.2;

%CU = CU*1000; % Compression for e format 20 oct 2005
%CV = CV*1000; % This multiplication is compensated by scale factor in source file
%fprintf(fp,'%s\n',['    ' sprintf('%10.5f',prd) '     ' sprintf('tsu%d',i) ' 1x1_mantle_k3.model ' sprintf('tsu%d',i) '.source' ' 1']);
fid = fopen('c:\manoj\ocean\agu07\ecco_annual_only_imag.source','wt');
fid1 = fopen('c:\manoj\ocean\x3dg\header.txt','rt');
%the following format for alexei's x3dg is tested ok 20 oct 05
           for ii = 1:34,
                s=fgetl(fid1);
                fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:180,
    for jn = 1:360,
       fprintf(fid,'%12.4e\n',CUr(in,jn));
   end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:180,
    for jn = 1:360,
        fprintf(fid,'%12.4e\n',CUi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
 for in = 1:180,
    for jn = 1:360,       
        fprintf(fid,'%12.4e\n',CVr(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
  
 for in = 1:180,
    for jn = 1:360,       
        fprintf(fid,'%12.4e\n',CVi(in,jn));
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:37,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
            fclose(fid);
            fclose(fid1);

fclose all;